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Stellar Turbulent Convection: A Selfconsistent Model 

V. M. Canute- \ I. Goldman 2 ' 3 - 4 and I. Mazzitelli 5 

ABSTRACT 

We present a selfconsistent model for stellar turbulent convection which is 
similar in spirit to the CM model (Canuto & Mazzitelli 1991) since it accounts 
for the full spectrum of the turbulent eddies rather than only one eddy, as done 
in the mixing length theory (MLT). The model differs from the CM model in 
the treatment of the rate of energy input n s (k) from the source that generates 
the turbulence. In the present model, n s (k) is controlled by both the source and 
the turbulence it ultimately generates, thus ensuring a selfconsistent modeling of 
the turbulence. This improves the CM model in which n s (k) was taken to be 
equal to the growth rate of the linear unstable convective modes. 

However, since the formulation of a selfconsistent treatment is far from 
simple, we were forced to use a representation of the nonlinear interactions 
less complete than the one in the CM model. The ensuing equations were 
solved numerically for a wide range of convective efficiencies. The results are 
the convective flux, the mean square turbulent velocity, the root mean squared 
turbulent pressure and the turbulent viscosity. 

We implemented the model in the ATON stellar structure code and computed 
the evolution of a solar model. The results are generally similar to those of the 
CM model and thus quite different from the MLT. The present model requires 
a smaller overshoot into the upper radiative zone than does the CM model, 
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in accord with recent empirical estimates. Application to POP II stars and 
comparison with the very metal-poor globular cluster M68 yields an age in the 
range 11 -j- 12 Gyr. This is somewhat younger than the CM age, which in turn 
is younger than the corresponding MLT age, a result of possible cosmological 
interest. 

Subject headings: convection — stars: interiors — stars: evolution — Sun: 
interior — turbulence 



1. Introduction 

Recently, Canuto & Mazzitelli (1991, CM) proposed an improved model for stellar 
convection. Being derived from a turbulence model, it takes into account the contribution 
of the full spectrum of the turbulent convective eddies, to the convective flux. In stellar 
interiors the microscopic viscosity is very small compared to the turbulent viscosity, 
implying that the turbulent spectrum spans many decades in wavenumber space. Therefore, 
in this respect, the CM model represents a significant improvement over the mixing length 
theory approach (MLT) which is a one eddy (the largest) approximation to the spectrum. 
Moreover, in the CM model, the turbulent mixing length scale is the depth z, so there is no 
need for an adjustable free parameter like the MLT a-parameter. The resulting convective 
fluxes are higher than those of the MLT for high convective efficiencies, and smaller than 
them for low efficiencies. The model performs better than the MLT when applied to stellar 
structure (D'Antona & Mazzitelli, 1994; Mazzitelli, D'Antona & Caloi 1995, Stothers 
& Chin 1995, Althaus & Benvenuto 1996), helioseismology (Baturin & Mironova 1995, 
Monteiro et al. 1995, Antia & Basu 1995) and stellar atmospheres (Kupka 1995). 

In the CM model, the turbulence spectral function E{k) is determined by the 
timescale controlling the energy input from buoyancy, that is, the source that generates the 
turbulence. This timescale is expected to depend on the parameters of the source as well as 
on the turbulence spectrum itself. The quantification of the latter dependence, within the 
CM model, is far from obvious. Thus, CM (1991) assumed that the above timescale can 
be approximated by the inverse of the growth rate of the unstable modes of the linearized 
equations. By construction, the latter depends only on the source and is independent of the 
turbulence it generates. The linear rate was used also by Canuto Goldman and Chasnov 
(1987, CGC) who, generalizing the work of Canuto and Goldman (1985), proposed a model 
for fully developed turbulence. The linear rate was also employed by Hartke Canuto & 
Dannevik (1988) in the framework of a DIA (Direct Interaction Approximation) model for 
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turbulent convection and by Canuto, Cheng, Hartke, & Schilling (1991) for EDQNM (Eddy 
Damped Quasi-Normal Markovian) models. 

The rate of energy input in the CM model must be improved since a fully developed 
turbulence is expected to regulate the energy input from source (buoyancy). The lack 
of feedback from the turbulence on the energy input, prevents the model from being 
selfconsistent. However, because of the complex structure of the EDQNM closure, the 
implementation of a selfconsistent rate into the formalism of the CM model is not simple. 
Thus, we were forced to simplify the structure of the nonlinear interactions, so as to be able 
to formulate a workable selfconsistent treatment. 

In modeling the nonlinear interactions, we follow the work of CGC (1987). However, 
in the present model we generalize the definition of the eddy correlation timescale, thus 
correcting some shortcomings in the physics involved, and leading to an improved closure. 
The resulting effective rate of energy input from the source (buoyancy) depends now on 
both the source and the turbulence. As stated, the present approach is complementary to 
that of the CM model. Here, the focus is on the selfconsistent rate of energy input and not 
on the closure, which is much simpler than that of the CM model. 

In spite of the highly nonlinear structure of the model equations, it is possible to 
solve them directly with no need for iterations. This reduced considerably the amount of 
numerical work required, and allowed for an exploration of the model for a wide range of 
values of the convective efficiency, S, defined in equation (51). For each value of S we 
obtained the spectral function which determines the turbulence bulk quantities. Thus, we 
computed, as functions of S, the convective flux, the turbulent viscosity, the turbulent 
mean squared velocity and the root mean squared turbulent pressure. 

We have applied the new model to the main sequence evolution of a solar model as well 
as to the evolutions of an extreme POP II chemical composition (Y=0.23, Z=10~ 4 ) stars 
with M < 0.9 M Q . The results are generally similar to those of the CM model. However, 
the new model has the advantage that the overshoot required to fit the solar model is much 
smaller, in accord with recent empirical estimates, and the ages of globular clusters are also 
somewhat smaller than the corresponding ages in the CM model (which in turn are smaller 
than those derived within the MLT framework). 



2. The Model 
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2.1. The Rate Controlling Energy Input from the Source 

Let us consider a fully developed stationary turbulent convection characterized by the 
spectral functions F(k) , G(k) and H(k), of the turbulent velocity, temperature fluctuations 
and the turbulent convective flux, respectively. Before doing so, we write the dynamical 
time-dependent equations obeyed by these spectral functions (Yamaguchi 1963). 

^F(k) + uk 2 F(k) = gaH(k) + T F (k) , (1) 
^G(k) + X k 2 G(k) = (3H(k) + T G (k) , (2) 

and 

® H (k) + (u + X )k 2 H(k) = (5r{k)F{k) + garG(k) + T H (k) , (3) 
at 

where T>, Tq, and Th denote the nonlinear transfer terms for the turbulent velocity, 
temperature, and convective flux. Here, g is the gravitational acceleration, a is the 
coefficient of thermal volume expansion at constant pressure (equaling T" 1 for an ideal gas), 
(3 is the superadiabatic temperature gradient, v is the microscopic viscosity, and x is t ne 
microscopic thermometric conductivity appearing in the expression for the conductive flux 

dT 

Fcond = ~ c vPX~^ ■ (4) 

In stellar interiors, the dominant conductive flux is the radiative flux and thus x is the 
radiative thermometric conductivity. The function r(k) is given by 

with x(k) = k\jk\ measuring the anisotropy of the eddy corresponding to the wavenumber 
k. Here, kh and k v stand for the horizontal and vertical wavenumbers, respectively (the 
vertical direction is that of the gravitational acceleration). 

In equation (1), gaH(k) plays the role of the energy source driving the velocity 
fluctuations. More precisely, it equals the rate of energy per unit mass and per unit 
wavenumber, fed to the turbulence velocity field at wavenumber k. The term vk 2 F(k) is 
the rate of energy per unit mass and unit wavenumber dissipated at k by the microscopic 
viscosity, while —T F (k) represent the rate of energy per unit mass and unit wavenumber 
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transferred to wavenumbers other than k. Analogous interpretations apply to equation (2) 
which describes the temperature fluctuations field. 

As stated above, we are interested in stationary turbulence. Thus, in equations (l)-(3) 
we set the time derivatives of the spectral functions equal to zero. Following CGC (1987), 
we assume that 



T H (k) TM/c) TgW 

H(k) F(k) G(k) 1 ' 

which equations (1), (2), and (3) show to be equivalent to the assumption: 

H{k)={r{k)F{k)G{k)) 1/2 . (7) 

For equation (7) to be satisfied, the velocity and temperature fluctuations, at any k, must 
be in phase. This is expected to be the case for turbulent convection where the temperature 
plays the role of an active scalar whose fluctuations drive the velocity fluctuations. 

We adopt the simplifying assumption that the transfer terms in equations (1) and (2) 
describe transfer from small to large wavenumbers only (from large spatial scales to small 
ones). While justified in three dimensional turbulence where energy flows predominantly 
from large to small scales, this simplifying assumption neglects non-local interactions (in k 
space) that give rise to reverse transfer (backscatter) . 

Denoting by k the wavenumber corresponding to the largest eddy, we shall represent 
the transfer as 



/* T F (k')dk' = -u t (k)y(k) (8) 



where 



y(k) = f F(k')k' 2 dk' (9) 

Jk 

is the mean squared turbulent vorticity and vt(k) is the turbulent viscosity at wavenumber 
k, exerted by all eddies of smaller size (larger k). As such, it is expressed as an integral 
with limits k and oo. We shall return later to the definition of the integrand. Similarly, 



where 



/ T G (k')dk' — —xt(k)w(k) 

Jko 



(10) 
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w (k) = f k G{k')k' 2 dk' (11) 

is the analog to the mean squared turbulent vorticity and Xt(k) is the turbulent conductivity 
at wavenumber k exerted by all eddies of smaller size (larger k). With these closures, we 
obtain by integrating equations (1) and (2) 

ga J*H(k')dk'=(v + v t (k)^y(k) , (12) 

and 

P £ H(k>)dk> = [x + X t{k)^w{k) . (13) 

The last two equations allow us to express the spectral function G(k) in terms of F{k). We 
obtain 

G(k) = —X(k)F(k) , (14) 
ga 

with 



A « = v m ' • (15 > 

E,(*) = V -±^§- , (16) 
X + Xt{k) 

where a prime denotes a differentiation with respect to k. Substituting equation (14) in 
equation (7) yields 



H{k)=[^-r{k)\{k)) ll2 F{k). (17) 



,ga' 

Thus, once F(k) is known so are G(k) and H(k), implying that one needs to solve only for 
the spectral function F(k). In order to do so, let us return to equation (12) and express 
H(k) in terms of F(k) using equation (17). We obtain 



J* (n s (k') + vk ,2 ^F{k')dk' = (v + u t (k)^y(k) (18) 
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where n s (k) is the shorthand abbreviation 



n s (k) = -vk 2 + [ga(3r(k)\(k)\ , (19) 

whose dimension are inverse time. The combination n s (k) + vk 2 plays, in equation (18), 
the role of the inverse of the timescale controlling the energy input into the turbulence at 
wavenumber k from the source (buoyancy in the present case). More specifically, the net 
rate of energy input from the source per unit mass and unit wavenumber, at k, is n s (k)F(k). 

It is obvious, from equation (19), that n s (k) depends on the source (buoyancy) and 
on the turbulent state. Thus, it conforms with the requirements discussed in §1. In what 
follows we show that equations (18) and (19) allow us to express n s (k) in terms of the 
turbulent viscosity u t (k). First, differentiate equation (18) with respect to k. The result is 

n.(k)-v t (k) , ^ = v t (k)k* . (20) 
Next, use equations (15), (20) and the definition of y(k), equation(9), to get 

m = Mk) + ^y{n s (k)-u t k 2 ). (21) 

With the help of equation (19), we transform equation (21) into a second order algebraic 
equation for \(k): 



X(k) = Z 2 t (k)4 + (ga(3r(k)X(k)) 1/2 k-^ , (22) 
whose positive solution is 




We adopt a relation between Xt(k) and u t {k) 

/ \ 1/2 

Xt(k) = (x 2 + ^ 2 v 2 {k)) -x (24) 

so that, as in CGC (1987), u t (k)/xt(k) equals the constant a t for large Pecle numbers, 
Pe = v t (k)/x, while u t (k)/xt(k) = 2a 2 x/^t(k), for small Pe. Using equation (24) to express 
S t and Xt in terms of v t , equation (23) now takes the form 
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with 



A(k) = (i+vr 2 v?(k) x - 2 ) 



3/2 



and 



= 1 



wt(k) 
°tX 2 



+ 



°tX 2 , 



+ A(k)- 



ga(3xe 2 T{k) 



2\ 1/2 



V 



With this X(k), n s (k) of equation (19) becomes 



(25) 



(26) 



(27) 



n s (fc) = — vk + 



2 , ga(5r{k)B{k) 



2 X k 2 A{k) 

Thus, we succeeded in expressing n s (k) in terms of the turbulent viscosity Vt{k), even 
though the latter is still unspecified. 



(28) 



2.2. The Eddy-Correlation Timescale 

In order to solve equation (18) (or equivalently equation (20)) we need an additional 
relation between the turbulent viscosity u t (k) and the spectral function F(k). Without loss 
of generality, v t {k) can be written as 

v t {k) = l°° ^rdk' (29) 
v ' Jk n*(k') v 1 

where n*(k) has the dimensions of an inverse time. In order to determine n*(k) let us focus 
on its physical meaning. Differentiation of equation (29) yields that the turbulent viscosity 
contributed by eddies in the wavenumber interval (k, k + dk) is 

du t = n* c (ky 1 F(k)dk . (30) 

Thus, n*(k) is proportional to the inverse of the eddy-correlation timescale at wavenumber 
k, or heuristically, the inverse of the timescale for the eddy breakup. The eddy is damped 
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because of two processes: interaction with the turbulent viscosity and the interaction with 
the source (microscopic viscosity is not considered here since in stars v << u t ). One can 
envisage the eddy as being "scattered" by "collisions" with the smaller eddies (turbulent 
viscosity) and by the source that drives energy into the eddy. The effective rate for the 
breakup (the inverse of the correlation timescale) will be taken to be the sum of the rates 
for the two processes as if they were operating independently: u t {k)k 2 + n s (k), plus the sum 
of the rates for each process which is now affected by the other. Between "collisions" due to 
one process there is a random walk due to "collisions" of the other process. Thus, these last 
two rates are n s (u t k 2 jn s ) l l 2 and v t k 2 (n s / \v t k 2 )) l l 2 . Summing up the four rates we obtain 



7 n*(A;) = v t {k)k 2 + n s {k) + 2^ t {k)k 2 n s {k)^ 12 = (iv t {k)k 2 ) 1 ' 2 + n^k) 1 ' 2 ^ (31) 

where 7 is a proportionality constant, determined by the normalization of y(k) in the 
inertial range. In this range equations (20), (29) and (31) yield 

y(k) = ^ (u t (k)k 2 ) 2 (32) 

while equation (18) results in 

e = y(k)v t (k) (33) 

with e = e(oo) where 

e{k) = J k k {n s (k') + isk' 2 ^jF(k')dk' = (v + v t (k))y(k) (34) 

is the energy rate per unit mass supplied to the turbulence from the driving source at all 
wavenumbers smaller than k. Combined together, equations (32) and (33) imply that 

y(k) = e 2/3 7 -l/ 3A; 4/3 (35) 

This should coincide with y(k) corresponding to the Kolmogorov spectral function 



Thus, we obtain 



</«(*) = \Ko<?^ . (36) 
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where K Q is the Kolmogorov constant. We used here the same normalization with respect 
to K as was used in CM (1991). 

With the definition of n*(k), equation (31), the model equations can be solved for F(k), 
and consequently for G(k) and H(k). We introduce for notational convenience (following 
the convention of CGC (1987)) the rate n c (k) defined by 

v t k 2 = 7 n c (/c) . (38) 

In the model for large scale turbulence of Canuto & Goldman (1985) n* was identified with 
n s , neglecting the role of the turbulence. In CGC (1987) n* = n c was used, neglecting the 
role of the source. Here, n* depends both on n s and on n c , as seen from equations (31) 
and (38), so it depends both on the source and on the turbulence. At k , n* = 4n c and 
therefore is 4 times larger than in CGC (1987). For high values of k (practically few times 
ko) n* — > n c . It is of interest to note that the inverse of the timescale for two-times velocity 
correlation according to the DIA model, indeed conforms to the present definition of n* (see 
figures 3 and 4 in Canuto & Battaglia 1988). In particular at k it is indeed ~ 4 n c . 

2.3. Differential Equation, Spectral Function, and Convective Flux 

From equations (29) and (38) we have 

F(k) = -^nl(k) , (39) 

and equations (20) and (38) lead to 

y(k) = n* c (k)(jn c (k)-n a (kj) . (40) 
Combining equations (39), (40) and (9) yields the differential equation for n c (k) 

2n* c (k)n' c (k) + n* c (k)'(n c (k) - 7 _1 ^(A;)) - T^n^KW ~ 2n c (k)n* c (k)^ = (41) 

with n s (k) and n*(k) defined as functionals of n c (k) through equations (28) and (31), 
respectively. A solution for n c (k) will also yield these two rates, as well as the spectral 
function F(k) and the mean squared turbulent vorticity y(k). 

The main objective of this work is the computation of the turbulent convective flux 
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roo 

F c = c p pu 3 9 = c p p / H(k)dk = c p pPx® • (42) 

Here 113 is the turbulent velocity in the flux direction, 9 is the temperature fluctuation, and 
the overline denotes ensemble averaging. The dimensionless convective flux $ is determined 
by equations (12) and (34) to be 

$ = eigaPxr 1 ■ (43) 

The value of e will be determined by using the fact that e(k) increases with k and reaches 
its asymptotic value already in the inertial range. Thus, we shall follow the solution of 
n c (k) and the corresponding e(k), from k up to to the inertial range, until e(k) saturates 
to e. From equations (34), (38) and (40) we find 

e(k) = n* c (k) [v + (W) - n.(*)) • (44) 

thus, once n c (k) is obtained, e(k) is determined too. 



3. Solution Procedure 

Turn now to the equation (41) which is a first order differential equation for n c (k) and 
thus a boundary value is required for a unique solution. Below we find the value of n c (k ). 
The value of k is determined by the width of the layer A, 

k = j(l + x y/ 2 (45) 

with xq = x(ko)- Since by definition F(ko) = and y(ko) = 0, equations (39), and (40) 
imply that 

jn c (k ) = n s (k ) = u t (k )ko (46) 

and 

(&-(£)!.">■ (47) 

The last two equations together with equation (28) determine n c (ko) 



-yn c (k ) = (Ttxklzo , 

where z is the solution of 



(48) 



<7 

h Zo 



(i + 4) 



1 + x l 



-cr 



-2 o2 



st 



(49) 



with 



ft = tT 4 (1 + xo) - ^ = 0.00456275 , (50) 
where S is the dimensionless product of the Rayleigh and Prandtl numbers 

S= 9 -^. (51) 

X 

S, or more precisely Si, is a measure of the convective efficiency. With the determination 
of n c (ko), equation (41) can be solved numerically to yield n c (k) for any value of S. 

Since we deal with stellar interiors for which a = ux^ 1 « 1, we take in all the 
equations the limit of v — > 0. Physically this means that the wavenumber at which the rate 
of dissipation by the microscopic viscosity becomes equal to jn^k), is much larger than fc 0) 
so that the spectrum exhibits an inertial range over many decades of k. For the limit of 
v — > 0, appropriate for stellar interiors, equation (49) yields an analytic solution for z as 
function of S, 



l + \/l + 4S?T 2 (lK 



implying that 

X 

in c \K G ) = n s {K ) = yyap) ■ r/ yD ) = 

with 



7 n c (A;o) = n s (k ) = (ga^^S) = ^S^S) (54) 



Vo{S) = I (55) 

/ i + v /i + 4s 1 V(ivr 2 

and r(l) = x (l + xq)' 1 . It is convenient to express the equations in terms of normalized 
rates 



Vs(k) = 



n s (k ) 



n s (k ) 



and 



7w;(fc) 

n s (k ) 

and to introduce a dimensionless wavenumber 



k 

q = T ■ 
k 

Using equations (26), (27) and (28) in equation (47) results in 

^r(k)\ 

We adopt, following CGC (1987), 



/ feo 



= 



A;A X 2 



x(k)=[— -1 



7T 



for which equation (60) yields 



k 

Thus 



3y/ 2 7r 

2J A 



and 



The differential equation for r) c (q) is obtained from equation (41) 
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WMrtM + V*MMq) - Vs(q)) ~ V*Mv' s (q) ~ 2vMv*M- q = , (65) 

where now the prime denotes differentiation with respect to q. By definition [see eqns. (56), 
(57), (58)] 



Vc (q = l) = Vs (q=l) = l ; v *(q = l)=4, (66) 
for any S. We have now 

Sl^TjqWq) 

^ = WW) (6?) 

where 
and 

S(?H1+|1 + 1^|MMV /2 . (69) 

The use of the normalized rates ensures that the computed quantities will not be very 
small (large) even for very small (large) S values, thus improving the numerical accuracy. 
Equation (44) becomes now, 

= ^(9aPx)vo s l /2 q' 2 Vc(q)v*(q) (vdo) - (70) 

the dimensionless convective flux is given by 

$(S) = ±rfiSl' 2 qJ 2 Vc(qM(q f ) (vdqf) - Vs(q f )) (71) 

where qf is the upper value of q for which equation (65) is solved, and is well inside the 
inertial range (thus, e(qj) = e). 
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4. Results 



4.1. Convective Flux 



We solved equation (65) for rj c (q), with the initial condition of equation (66), for S in 
the range 1CT 4 — 10 20 . We adopted the commonly used value of a t = 0.72. Each solution 
yields r] s (q) and r]*(q) and the spectral function F(q). For each such solution we computed 
e(g), followed it to saturation, and thus obtained e(S) and $(5*). 

In Figure 1 we show e(q) in units of (K /l.5) 3 ga(3xi f° r S — 10 6 . The qualitative 
behavior of e(q) is typical for any value of S: it starts from zero and saturates in the inertial 
range to e. From equations (70) and (71) it follows that the asymptotic value of the graph 
equals <£> in units of (K /1.5) 3 . 

In Table 1 we list ®(S) (rounded to 4 figures) for 20 representative values of S. From 
equation (71) it follows that <&(S) oc 7 _1 oc Kq. The values shown in Table 1 are in units of 
{K /l-5 f. The limiting behavior of $(S) is given by 



2.65 x 10~ 5 ( ^ 
V 1.5 



S 2 



S« 1 , 



(72) 



and 



In applications of $(S) to stellar structure codes, it is useful to have an analytic fit 
formula to the convective flux. We derived such a fit with a deviation ^ 3%: 



^ = F 1 (S)F 2 {S) (74) 

where 



F 1 (5) = (^) 3 a5*((l + 65) m -l) B (75) 



'K f 

*W) = ( 
where 

a = 10.8654 , b = 0.00489073 , k = 0.149888 , m = 0.189238 , n = 1.85011 



and 
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Table 1. The Results of the Model 



c 




$ 




9 

V 


Pt 




10~ 4 


2.65 x 10~ 13 


0.3091 


2.814 


2.112 x 10~ n 


5.3408 x 10~ 12 


0.2529 


0.01 


2.65 x 10~ 9 


0.3092 


2.814 


2.112 x 10~ 7 


5.341 x 10~ 8 


0.2529 


1. 


2.65 x 10~ 5 


0.3148 


2.819 


2.112 x 10~ 3 


5.341 x 10~ 4 


0.2529 


10 


2.649 x 10~ 3 


0.3663 


2.863 


0.2111 


0.5339 


0.2529 


30 


2.38 x 10~ 2 


0.4811 


2.96 


1.895 


0.4793 


0.2529 


10 2 


0.2557 


0.8745 


3.203 


20.52 


5.186 


0.2527 


300 


1.88 


1.806 


3.409 


157.6 


39.63 


0.2515 


10 3 


10.23 


3.334 


3.045 


1.019 x 10 3 


252.4 


0.2478 


10 4 


85.01 


5.82 


2.065 


1.685 x 10 4 


4.017 x 10 3 


0.2384 


10 5 


388.5 


7.381 


1.59 


1.954 x 10 5 


4.515 x 10 4 


0.2311 


10 6 


1.442 x 10 3 


8.312 


1.348 


2.049 x 10 6 


4.648 x 10 5 


0.2268 


10 7 


4.917 x 10 3 


8.849 


1.209 


2.085 x 10 7 


4.680 x 10 6 


0.2245 


10 8 


1.615 x 10 4 


9.155 


1.125 


2.094 x 10 8 


4.687 x 10 7 


0.2238 


10 9 


5.21 x 10 4 


9.326 


1.072 


2.099 x 10 9 


4.688 x 10 8 


0.2233 


10 10 


1.665 x 10 5 


9.42 


1.038 


2.101 x 10 10 


4.688 x 10 9 


0.2231 


10 12 


1.679 x 10 6 


9.498 


1.001 


2.101 x 10 12 


4.688 x 10 11 


0.2231 


10 14 


1.684 x 10 7 


9.528 


0.9863 


2.101 x 10 14 


4.689 x 10 13 


0.2231 


10 16 


1.685 x 10 8 


9.534 


0.9795 


2.101 x 10 16 


4.689 x 10 15 


0.2231 


10 18 


1.685 x 10 9 


9.534 


0.9764 


2.101 x 10 18 


4.689 x 10 17 


0.2231 


10 20 


1.685 x 10 10 


9.534 


0.9751 


2.101 x 10 20 


4.689 x 10 19 


0.2231 



$ is in units of (Ko/l-5) 3 and its high S limit scales as (o-f/0.72) 3 / 2 . v 2 and pt are in units 

of (|) 2 (ff) 3 and p(t9 2 (t5) 3 > respectively. Their high S limits scale as (<r t /0.72). The 
dimensionless ratio, C, is defined in equation (88). 
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F 2 (S) = 1 + 



1 + dS" 



1 + fff 



(76) 



where 



c = 0.0108071 , d = 0.00301208 , e = 0.000334441 , / = 0.000125 , 



p = 0.72 , q = 0.92 , r = 1.2 , t = 1.5 . 



In order to judge the quality of the fit, in Figure 2 we display the ratio between the fit 
function and an interpolation of the numerical values, as a function of log S. 

Also shown in Table 1 are the ratios $(S)/$mlt(S) and ®(S)/$cm(S) where $mlt 
and $cm are the values corresponding to the MLT, 



and to the CM model (their eq. [32]). Figures 3 and 4 display &(S)/§mlt(S) and 
&(S)/$cm(S), respectively. Note that the qualitative behavior of &(S)/$mlt(S) is similar 
to that of the CM model — higher flux for high S and lower flux for low S values. The 
comparison with $cm(S) shows that while the two models yield essentially the same flux 
for high S values, the new model predicts higher fluxes for intermediate and low S values, 
and the flux ratio is maximal for S ~ 300. Comparisons of $(S) computed within the new 
model with $(S) computed in a model with the same definition of n* but with n s equal 
to the linear growth rate, indicates that the above local maximum at S ~ 300 is a feature 
resulting from the use of the selfconsistent rate. We recall that Si = 0.00456275* rather 
than S is the measure of the convective efficiencies. Thus, the borderline between low and 
high efficiencies is around S = 300, which is also where the ratio of the new convective flux 
to that of the CM model is maximal. 



4.2. Turbulent Velocity, Turbulent Pressure and Turbulent Viscosity 

The mean squared turbulent velocity is defined by 




(77) 




(78) 



which with the use of equation (39) can be expressed as 
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1 t°° 

v 2 = X 2 A- 2 l- 1 ^ r r—-r Sv 2 o(S) f(q)dq, (79) 

with 

/(?) = -»?:(?) (^)' • (so) 

We computed v 2 for each value of S and the results are presented in Table 1. The limiting 
behavior is given by 

3 /,,\ 2 



v 2 = 2.1117 x 1(T 3 (^) (£) ; 5«1 (81) 



and 



" 2 = 21014 HifJ U#AxJ s ; s>>1 < 82 > 

The turbulent pressure is of importance in helioseismological models. Batchelor (1953) 
derived an expression for the mean squared turbulent pressure (for the case of isotropic 
turbulence) in terms of the spectral function F{k) 



p 2 = -p 2 ^ ^ F{k)F{k')I{k/k')dk'dk (83) 
where p is the mean density and the dimensionless integral I(x) is 



I(x) = I(l/x) = \(x 2 + x- 2 ) -\-\{x + x- 1 ) (x - x- 1 ) 2 lnl±^_ . (84) 

I(x) — > for x — ■> , oo and is maximal at x — 1, where it equals 2/3. Thus, the pressure is 
mostly contributed by k! ~ k which are close to the maximum of F(k). Using equation (39) 
we obtain 

Pt = ^fP 2 X^- \ i{1 + Xo)2 S 2 V 4 o(S) I I f(q)f(q')I(k/k')d q 'd q . (85) 

The computed values of the root mean squared turbulent pressure are displayed in 
Table 1. The asymptotic behavior is given by 
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p t = 5.3408 x lO-^pf^) 2 (^) 3 S 2 ; S«l (86) 



AJ V 1.5 
and 

*=°-^(x)' ; s » i - « 87 > 

A quantity of interest is C(S) defined as 

C(S) = ^= (88) 

The computed values of C are listed in Table 1. As can be seen, it is almost constant 
for all values of S, ranging from 0.253 for low S to 0.223 for high values of S. 

As with the convective flux it is useful to have analytical fit formulae for the mean 
squared turbulent velocity and for the root mean squared turbulent pressure. We derived 
such fits which represent the numerical values with precision better than 3%. For the mean 
squared velocity we derive the fit 

1,2 = (I)* F3iS) FiiS) • m 

,Kq}> 0.00101392^ 
V1.57 1 + Vl + 0.000017848 S 2 

and 

p (o\ — a oo 8 oo , 2.256815 (-1. + 0.000777055 S -"*™) 

F 4 (S) - 6.39899 + i. + .000777055 ■ ( 91 ) 

Similarly for the root mean square of the turbulent pressure we find 

Pt = p(j) 2 F 3 (S)F 5 (S) (92) 



where 



with F 3 (S) given by equation (90) and 



F. {S) = 1,^8 + OA^ -^—p^ ■ m 
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Finally, the turbulent viscosity, is given in the present model already in an analytic 
form. From equations (38), (44), (45), (54) and (55) we find that 

„ , 0.00215086 5 , . 

v t = u t (k ) = \ -, = . (94) 

VI + VI + 0.000017848 S 2 

The limiting behavior of v t is given by 

u t = 0.00152089x5 ; S « 1 (95) 

and 

v t = 0.03309 x S 1/2 ; S » 1. (96) 



5. Application to stellar models 

The new convective fluxes have been included in the ATON stellar structure code (for 
an update on the physical and numerical details of the code, see Mazzitelli et al. 1995 and 
references therein). We have computed the main sequence evolution of a solar model as well 
as a set of evolutions for Pop II stars having M < O.9M , from the zero age main sequence 
to the base of the red giant phase. 

Before turning to a detailed discussion of the results, we recall that, the turbulent 
length-scale A at a given depth z inside a convective region, must also include the thickness 
OV of the overshooting layers (if any) beyond the formal Schwarzschild boundary (see 
DAntona and Mazzitelli 1994). At present, the OV phenomenon has not yet been fully 
quantified in a reliable way (Umezu 1995) even though the underlying equations have been 
derived (Canuto 1993). However, empirical evidence from comparisons between stellar 
models constructed with local convection theories, and observations of intermediate mass 
main sequence stars in young open clusters, suggests quite stringent limits on the extent of 
the OV, namely, < OV < 0.2H P (Stothers & Chin 1992). Lacking a formal theory, we 
shall write 



A = z + a*Hl op . (97) 

where H t ° p is the pressure scale height at the upper boundary of the convective layer 
determined by the Schwarzschild criterion and a*, which should not exceed 0.2, can be 
regarded as a fine-tuning parameter. 
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We stress that the role of a* in the CM and in the present models is radically different 
from the role of the parameter a in the MLT model where A = aH p (z). In fact, a is a free 
adjustable parameter through which modelists try to capture all the physical uncertainties 
(e. g. opacities, convection, thermodynamics etc.) relevant to the evaluation of the effective 
temperatures of stars . Unfortunately, not only does this procedure seriously hinder the 
predictive power of stellar modeling (a is determined a posteriori), but the fit to the 
observed surface temperature does not automatically guarantee that also the internal 
temperature profile is correct. In this context, it is worth noting that Baturin & Mironova 
(1995) and Monteiro et al. (1995) have recently shown that the solar internal temperature 
profiles predicted by the CM model are in better agreement with helio-seismological data 
than those derived from the MLT. Moreover, Gabriel (1995) has shown that the MLT 
could be made to predict a CM-like internal temperature profile, provided that a is forced 
to vary inside the convective region, in a manner that represents an a posteriori fitting. 
This quite clearly shows that a) the MLT has no predictive power and b) the degree of 
artificiality that is required to make the MLT yield results that the CM model produces 
quite naturally. The parameter a*, on the other hand, quantifies a well identified physical 
process, the convective overshooting OV and, as seen in the following, in the present model, 
only a marginal amount of tuning is allowed anyhow. 

Finally, the convective flux in Table 1 is normalized to a value of the Kolmogorov 
constant of K Q = 1.5, as in the CM model. Since recent experimental data suggest higher 
values of K Q (up to ~ 1.9), we employed for the stellar modeling a fiducial value of 
K = 1.7 in the fit-formula, equation (75). Since, as discussed in §4.1, the convective flux 
scales as Kq, the flux used is a factor of (1.7/1.5) 3 larger than the numerical values in Table 
1. In Mazzitelli et al. (1995) K Q = 1.8 was used. 

Using the ATON code, and updating the low-T opacities according to Alexander and 
Ferguson (1994), we obtained a fit to the observed solar radius and luminosity at an age 
~ 4.55 Gyr (Bahcall 1989), and with a metal abundance Z=0.0175 (Grevesse and Noels 
1993), with Y~0.27 and a* ~ 0.08. The latter corresponds to a very small amount of 
overshooting of a few kilometers, and is in full accordance with the observational limits of 
Stothers & Chin (1992). On the other hand, had we employed the original CM model the 
required value would have been a* ~ 0.2 which is a borderline value. To find the maximal 
variance in T e ff allowed by the a* parameter, we computed two solar models with the 
borderline values a* = and a* = 0.2. The difference in T e ff between these last two models 
turned out to be <4%. 

Figure 5 shows the internal profiles of the dimensionless temperature gradient, 
d\ogT/d\ogP, in the region of the overadiabaticity peak for solar models computed with 
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the present, the CM and MLT models (for the MLT model a = 1.55). The results of the 
present model are very similar to those of the CM model but are quite different from the 
MLT results. The similarity of the results of the present and CM models, in spite of the 
difference (factor ~ 3) in the value of the convective flux for intermediate and low convective 
efficiencies, can be readily understood. In the more external convective layers, the density 
is quite low (p < 1CT 7 ), and LogS* < 2. The turbulent flux is also quite low, and convective 
energy transfer is very inefficient. Most of the flux in this region is therefore carried by 
radiation, and the temperature gradient sticks to the radiative one. In deeper convective 
layers, where p, S and A are larger, convection begins to be efficient and the value of the 
temperature gradient is determined by the turbulent convective flux. However, for S values 
such that Log S ^> 3, the new fluxes are very close to the CM ones, and so is the resulting 
temperature gradient. The vicinity of the gradient peak, LogS* ~ 2, is close to the S value 
where the present and CM models differ the most ( see Fig. 3). Thus, this is the region 
where we can expect some sizable difference between the temperature profiles, as Figure 5 
indeed shows. 

Because of the similarity of the results from the present and the CM models, the 
experimental benchmark provided by helioseismological data cannot discriminate between 
the two models (Antia and Basu 1995). However, on the basis of stellar modeling, we stress 
that the new fluxes require a lower value of a* to fit the sun, which is more in agreement 
with the results of Stothers & Chin. Since the low-T radiative opacities are probably still 
slightly underestimated, and larger opacities require a larger a* to fit the sun, we prefer 
the present new fluxes over the original CM values. The reason is that the latter, once the 
updated values of low-T opacities become available, could require values of a* larger than 
allowed by the observational upper limit on overshooting. 

As a further check, we have also applied the new fluxes to the computation of 
evolutionary tracks and isochrones for stars with Y = 0.23 and Z = 10~ 4 . The isochrones 
are shown in Figure 6, together with the fiducial sequence of the globular cluster M68. 
Details on the computations of both tracks and isochrones, on the observational to 
theoretical correlations, as well as the chemistry, reddening and distance modulus for M68 
can be found in Mazzitelli et al. (1995). Here we simply recall that the apparent "kinks" in 
the isochrones are a true physical feature, which is expected and explained within the CM 
framework, and which exists also in the present model. 

The age of the cluster is in the range 11 -j- 12 Gyr, somewhat younger than the 12 13 
Gyr found with the CM fluxes, which itself is younger than the 13 -j- 15 Gyr derived within 
the MLT. Whether this difference is significant for solving the age conflict of globular 
clusters with the age of the universe, following from the recent determinations of high values 
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of H (Freedman et al. 1994, Pierce et al. 1994), is beyond the purposes of the present 
paper. This point will be discussed elsewhere (Canuto et al. 1996). 

In conclusion, the application of the convective fluxes of the present model to stellar 
structure does not alter appreciably the main results of the CM model. Nevertheless, the 
new model is preferable on two grounds. From the theoretical viewpoint, it determines the 
rate of input of energy from buoyancy to the turbulence in a more physically consistent 
manner. From the astrophysical point of view, it requires a smaller extent of overshooting, 
which is in better agreement with recent observational results. 

6. Discussion 

We have presented a selfconsistent model for turbulent convection based on a simplified 
treatment of the non-linear interactions among the eddies. The important novel feature 
of the present model is the formulation of a selfconsistent rate for energy input from the 
source (buoyancy) into the turbulence, which depends both on the source parameters and 
on the turbulence itself. This represents an improvement compared to the CM model where 
the rate of energy input was the growth rate of the linear unstable modes. The focus of the 
present model is on the selfconsistent rate of energy input at the expense of a less complete 
treatment of the non-linear eddy interactions. The latter is much simpler than in the CM 
model and describes transfer only from small to large wavenumbers. This representation 
neglects non-local (in k space) interactions between the eddies that lead to a reverse transfer 
(backscatter) that are included in the CM model. 

We have explored the model for a wide range of convective efficiencies and computed, 
numerically, the dimensionless convective flux, the turbulent squared velocity and the 
root mean squared turbulent pressure, as functions of the convective efficiency S. The 
results were fitted by analytical formulae with precision better than 3% over the range 
S = 10~ 4 -J- 10 20 . The turbulent viscosity in the model is already given by an analytic 
expression. The convective flux, is larger than that of the MLT for high convective 
efficiencies and lower than it for low convective efficiencies. This general behavior is similar 
to that of the CM model. The high S fluxes are very close to those predicted by the CM 
model but the intermediate and lower S fluxes are larger than those of the CM model. 

It is of interest to note that even the very simple CGC (1987) model shares the same 
qualitative behaviour of the convective flux relative to that of the MLT flux, as function 
of the convective efficiency The fact that three models differing in their treatment of 
the energy input rate and in the modeling of the nonlinear transfer, still yield the same 
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behaviour is quite intriguing. It suggests that accounting for the full turbulence spectrum 
(common to all three models) is far more important than the detailed way in which thye 
latter is done. 

We have applied the new model to the main sequence evolution of a solar model as 
well as to evolutions of POP II stars with M < 0.9 M . The convective turbulent length 
scale was taken equal to the depth in the convective zone, as in the CM model. The 
results are generally similar to those of the CM model. However, the new model has the 
advantage that the overshooting required to fit the solar model is smaller than in the CM 
model. Also, the ages of globular clusters are smaller than the corresponding ages in the 
CM model by ~ 1 Gyr, which may help alleviate a possible conflict between the ages of 
globular clusters and a high value of H . As already noted in §5, the similarity between 
the temperature profiles for the solar model, predicted by the present and the CM models, 
renders the models practically indistinguishable by helioseismological data (Antia & Basu 
1995). The situation is similar with regard to solar atmosphere modeling (Kupka 1995). 
However, atmospheres of cool stars are expected to yield observable differences between the 
two models (Kupka 1995). 

From the theoretical perspective, a more complete model is one that incorporates a 
selfconsistent rate of energy input while keeping the non-linear interactions in their full 
generality, as done within the CM model. Work in this direction is in progress. 

Finally, the present and the CM models are based on two-point correlations of the 
turbulent quantities. This methodology, preferred by the physics community, yields 
information about the spectral properties of the turbulence. Yet, its applicability to 
inhomogeneous and anisotropic cases is limited. An alternative approach, based on 
one-point correlation functions, is widely used in the engineering community and can handle 
anisotropy and inhomogeneity. While the spectral information is lost in this Reynolds Stress 
formalism, the method is easy to apply to non-local and space-dependent problems and 
has the potential to treat stellar convective overshooting (Canuto 1993). The method has 
already been successfully applied to the planetary boundary layer (PBL) which is the seat 
of strong convection (Canuto et al. 1994a) as well as to study the interaction between shear, 
vorticity and buoyancy at the surface of the sun (Canuto et al. 1994b). Work is in progress 
(Gabriel 1996, Houdeck 1996) to apply the same method to the study of helioseismology. 

This work was supported by the US-Israel BSF grant 94-00314, to I. Goldman 
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Fig. 1. — e(g) in units of (K /1.5) 3 gaj3x> f° r S — 10 6 . The asymptotic value shown is 
actually $ in units of (Kq/1-5) 3 . 
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Fig. 2. — The ratio between the fit function for 3>(S), eq. (74), and an interpolation of the 
numerical values of &(S). 
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Fig. 4.— The ratio $(S)/$ C m(S). 
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Fig. 5. — The dimensionless temperature gradient, d log T/d log P, versus pressure in the 
upper convective layer of the sun. The solid line corresponds to the present model. The 
dotted line corresponds to the CM model. The MLT (with a = 1.55) yields quite different 
results represented by the dashed line. 
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Fig. 6. — Isochrones in the HR diagram computed with the present model for an extreme 
Pop II chemical composition (Y = 0.23 and Z = 10~ 4 ). The squares mark the fiducial Turn 
Off region for the very metal-poor globular cluster M68. 



